Global transmission suitability maps for dengue virus transmitted by Aedes aegypti from 1981 to 2019

Mosquito-borne viruses increasingly threaten human populations due to accelerating changes in climate, human and mosquito migration, and land use practices. Over the last three decades, the global distribution of dengue has rapidly expanded, causing detrimental health and economic problems in many areas of the world. To develop effective disease control measures and plan for future epidemics, there is an urgent need to map the current and future transmission potential of dengue across both endemic and emerging areas. Expanding and applying Index P, a previously developed mosquito-borne viral suitability measure, we map the global climate-driven transmission potential of dengue virus transmitted by Aedes aegypti mosquitoes from 1981 to 2019. This database of dengue transmission suitability maps and an R package for Index P estimations are offered to the public health community as resources towards the identification of past, current and future transmission hotspots. These resources and the studies they facilitate can contribute to the planning of disease control and prevention strategies, especially in areas where surveillance is unreliable or non-existent.

The spatiotemporal dynamics of MBVs are primarily governed by the interplay among three factors: the physiological interactions between virus and vector, that between virus and host, and the population dynamics of the vector and host. The carrying capacity and seasonal oscillations of mosquito populations are influenced by many environmental and ecological factors including temperature [14][15][16] , humidity 16 , host population density 16,17 , vegetation types 18 and the presence of water sources 16,17 . Climatic factors are particularly important because they also alter each mosquito's potential to transmit the virus to new hosts by causing changes in physical and behavioral traits such as life span 19 , incubation period 20 and biting rate 21 . With data from experimental studies, it is possible to mathematically characterize these relationships between climatic variables and mosquito-viral traits and incorporate them into mechanistic transmission models 22 .
The basic reproductive number (R 0 ), which measures the average number of secondary infections generated by a single infectious host in a fully susceptible host population, is an important measure of transmission potential. For MBVs, expressions of R 0 involve many interacting host, pathogen and vector variables, some of which are difficult to parameterize due to limited data (e.g. geographical mosquito distribution and density). To overcome this challenge, different measures of transmission potential, called suitability measures, have been developed 18,21,23,24 . Vectorial capacity, one of the most widely known suitability measures, uses purely entomological variables to estimate the potential number of infectious bites that would arise from a single infectious person over one day 23 . A number of studies have built upon the concepts in vectorial capacity to develop new measures of transmission potential that consider the interactions between climate and mosquito-viral traits 22,24,25 .
Here, we update and apply a previously developed mosquito-borne viral suitability measure 24 , referred to as Index P, to estimate geographical maps and time series for climate-driven DENV transmission potential of the Ae. aegypti vector. Index P is a proxy for the transmission potential of a single adult female mosquito under conditions where susceptible hosts, the virus and its vectors are assumed to be present 26 . Index P is unlike other previously developed environmental suitability measures 18,25 in terms of its unique biological interpretation. It also accounts for additional host-pathogen factors that are included in classic transmission modelling frameworks (e.g. human infectious period), offering ways of directly parameterizing transmission models. It also builds upon purely temperature-dependent suitability measures 22 by incorporating the effects of humidity on entomological factors such as mosquito mortality rate. Local temperature and humidity time series are its main inputs, making its general framework unique in that it can be readily applied to any location for which such climate data exists and for any MBV for which there is empirical data on the relationship between climate and vector-viral traits. It has been successfully used to characterize the transmission potential and epidemiology of West Nile virus in Israel 27 , Portugal 28 , Brazil 29 , and Italy 30 , CHIKV and ZIKV in the Dominican Republic 5 and Mexico 31 , and DENV in Myanmar 26 , Brazil 24,32 and Mexico 31 . So far, the application of Index P has been performed at a small scale driven by specific and limited research goals. Given the sparsity of detailed epidemiological and ecological data for DENV (e.g. infection incidence and mosquito density) and the disproportionate impact of DENV in developing countries, there is a need to extend these analyses globally. Hence, we offer a database of spatiotemporal maps of Index P for 186 different countries and territories and an easy-to-use R package for Index P estimation as ready-to-use resources for the visualization and analysis of climate-driven DENV transmission potential globally over the last four decades.
In combination with previous work to map the distribution of DENV vectors 18,33,34 , these maps improve our understanding of past and current DENV transmission potential and have the capacity to inform disease control and prevention strategies. For example, in Central Africa where Ae. aegypti is predicted to be widely disseminated 18 but reported DENV cases remain low 35 , these maps can help determine whether the apparent mismatch is due to low transmission potential of the vector or under-reporting and poor surveillance. This information can then be used to identify high-risk, low-surveillance areas where seroprevalence surveys should be targeted. It thus addresses a central limitation of occurrence-based predictions of dengue infection risk 8 , which can underestimate risk in areas where reported incidence is low or absent due to limited surveillance rather than low transmission. The maps may also be used to identify regions where Ae. aegypti might not be present but due to high DENV transmission potential there is a high likelihood of spillover of sylvatic DENV into human populations 36,37 .

Methods
Fundamental theory of Index P. Index P is a climate-driven suitability measure for mosquito-borne viruses derived from a mechanistic model previously developed to study the transmission dynamics of Zika 38 and dengue 26 . The equation for R 0 in this model (Eq. (1)) can be decomposed into two components: the number of female mosquitoes per host (M) and the transmission potential of each female mosquito (P). Index P measures transmission potential solely based on P (Eq. (2)), which is interpreted as the reproductive potential of a single adult female mosquito in a fully susceptible host population. A description of the parameters in the equation for Index P can be found in Table 1. Given the availability of empirical studies that quantify the relationships between meteorological variables and viral-vector traits (e.g. extrinsic incubation period), it is possible to define equations for several parameters in the expression for P in terms of temperature (t) and humidity (u). Temperature and humidity time series data can be input into these equations to derive time series of the entomological parameters, which when combined with prior information on the selected vector/host/virus system can provide estimates of Index P over time. We have made two major updates to the original methodology of Index P 24 based on a revision of available empirical data sources. These updates include introducing a temperature-dependent function for the probability of transmission from human to mosquito per bite and a new temperature-dependent probability distribution of the extrinsic incubation period. ( ) Proposal of informative distributions of the biological parameters. The estimation of Index P requires some prior knowledge of the biology of the selected host, vector and viral system. In this study, we defined probability distributions for human hosts, Ae. aegypti vectors and dengue virus. Several mosquito species in the genus Aedes are known to transmit DENV to humans. We focused on Ae. aegypti because it accounts for the majority of vector-human transmission 39 , and it is the species for which the most empirical data exists on the relationship between climate and vector-viral traits. After a review of relevant literature, we defined probability distributions for six host-virus and vector-virus parameters: extrinsic Ae. aegypti-DENV incubation period, adult Ae. aegypti lifespan, adult Ae. aegypti biting rate, human lifespan, intrinsic human-DENV incubation period, and human-DENV infectious period (Table 1). These probability distributions are either directly sampled (μ h , γ h , σ h and t v ( ) γ ) or used as likelihood functions for the estimation of climate-dependent probability distributions and a u v ( ) ).

Climate-dependent functions for entomological parameters.
Weather-dependent functions are defined for each of the entomological parameters in the expression for Index P (Eq. (2)). The adult human mortality rate (μ h ), intrinsic incubation period (1/γ h ) and the human infectious period (1/σ h ) are taken to be climate-independent. As detailed in Obolski et al. 24 , adult vector mortality (μ v ) and the probability of transmission from vector to human per bite (φ v→h ) are modeled by temperature-dependent functions estimated in experiments on laboratory strains of Ae. aegypti by Yang et al. 40 (Eq. (3)) and Lambrechts et al. 41 (Eq. (4)), respectively. In the original implementation of Index P, the probability of transmission from human to mosquito per bite (φ h→v ) was assumed to be constant. To account for the temperature dependence of φ h→v and to capture threshold effects at extreme temperatures, we expanded the framework by modeling φ h→v with a ramp model estimated by Lambrechts et al. 41

Parameter Symbol Distribution Mean, SD Units Sources
Adult Ae.  www.nature.com/scientificdata www.nature.com/scientificdata/ Following the original approach by Obolski et al. 24 , the influence of humidity on adult mosquito mortality rate (μ v ) and mosquito biting rate (a v ) is introduced according to the following expressions. The time series for relative humidity is normalized to [0,1] with the ecological variables centered around the local average (u).
The total effect of climate on the entomological parameters is then modeled by combining the estimated temperature-dependent and humidity-dependent functions.
The scaling coefficients η and ρ in the expressions for mosquito biting rate (Eq. (11)) and adult mosquito mortality rate (Eq. (8)) are estimated from the user-defined probability distributions of a v and μ v , respectively. These distributions permit deviations from the ideal laboratory conditions used in the original empirical studies. The multiplicative coefficient η determines the magnitude by which the effect of temperature on the adult mosquito mortality deviates from that observed under laboratory conditions. The effect of temperature under field conditions is identical to that under laboratory conditions when η = 0, stronger when η > 1 and weaker when η < 1. The exponential coefficient ρ modulates the relative influence that humidity has on adult mosquito mortality and mosquito biting rate: humidity has no effect when ρ = 0 and a stronger effect when ρ > 0. For further information on the biological interpretation of these coefficients and their estimation, please refer to Obolski et al. 24 .
Using recent, high-resolution empirical data, we also updated the probability distribution of the extrinsic incubation period (EIP) of Ae. aegypti mosquitoes. We introduced a temperature-dependent probability distribution of γ v derived from observational studies of EIP in human subjects experimentally infected by wild Ae. aegypti mosquitoes 20 . Chan and Johansson 20 estimate the temperature-dependent distribution of the EIP of Ae. aegypti mosquitoes (1/γ v ) by directly fitting a log-normal time-to-event model to EIP observational data: Where α(t) and β(t) are temperature-dependent functions of the mean and standard derivation of the natural logarithm of the EIP ( Table 1).
Estimation of Index P. Index P estimation is based on a Markov chain Monte Carlo (MCMC) sampling framework which derives the distribution of the Index P time series using the climate-driven functions for the entomological parameters (μ v , φ v→h , φ h→v , a v ), the user-defined biological probability distributions (Table 1) Estimation and sampling of Index P globally. We used monthly temperature and humidity data, as well as informed probability distributions of the biological parameters as inputs to estimate the monthly Index P time series for each spatial pixel (~28 km 2 ) of 186 individual countries and territories from 1981 to 2019. For each country, spatial polygons were used to curate the average monthly surface air temperature and surface air relative humidity data sourced from Copernicus.eu 42 for all spatial pixels within its geographical boundary. The spatial pixel resolution used was 0.25° × 0.25° (~28 km 2 ), the maximum resolution available for the Copernicus.eu dataset.
The following steps were then performed for each spatial pixel. Using the climate data and the user-defined probability distributions for mosquito biting rate (a v ) and adult mosquito mortality rate (μ v ), we estimated the (2023)  www.nature.com/scientificdata www.nature.com/scientificdata/ climate-dependent posterior distributions of these two parameters by fitting the scaling coefficients η and ρ. For each MCMC chain, we drew 20,000 samples with the first 20% of samples treated as burn-in. Given a lack of information on the possible estimates of ρ and η, we assumed uninformative, flat priors for the two scaling coefficients in the ranges 0-10 and 0-20 respectively. We then drew 1,000 samples from the climate-independent distributions of μ h , γ h and σ h , 1,000 samples from the climate-dependent distributions of the γ v , μ v and a v time series, and solved the deterministic temperature-dependent equations of φ v→h and φ h→v . These sampled values were plugged into the equation for Index P, resulting in 1,000 monthly Index P time series. Finally, we calculated various summary statistics of the distribution of monthly Index P time series for each spatial pixel. Details on the calculated summary statistics are provided in the Data Records section. These estimates were collated at the country/territory level into TIFF files 43 .

Data records
For each country or territory, we calculated epidemiologically relevant summary statistics for the spatiotemporal distribution of Index P. These include: spatiotemporal map(s) for each month (468 layers), each year (39 layers), a typical year (12 layers) and the entire period (1 layer); spatiotemporal maps of the number of months Index P is above 1.0 for each year (39 layers), and for a typical year (12 layers); and spatiotemporal maps of the peak and trough timing during a typical year (12 layers each). The Index P spatiotemporal maps are projected with World Geodetic System (WGS)84, latitude/longitude coordinate system (EPSG: 4326; https://epsg.io/4326). The complete dataset and a description of its contents are available from a figshare repository 43 . The complete dataset is offered as a single Zip file which contains individual folders for each of the 186 countries or territories. Each folder contains a collection of spatiotemporal maps in the TIFF file format. The files are named according to the country/territory name (e.g. Brazil), the temporal resolution (e.g. monthly), and the relevant summary statistic (e.g. mean). Full details on the file names are provided in a separate content description file 43 . The Index P maps can be processed with packages such as raster 44 in the R programming language. A sample of the Index P maps and time series is in Fig. 2.  46 , Puerto Rico (1987-2012) 46 and Thailand (2003-2019) 47 , we performed validation analyses of Index P as a climate-based measure of transmission potential. We focused on these countries for several reasons: (1) DENV is endemic with high disease burden 8,48 , (2) the countries span diverse climate conditions with both temperature and tropical areas 49 , (3) there is large variation in the intensity and seasonality of observed DENV dynamics within Brazil 50 and (4) high resolution DENV case series are available. In accordance with the available spatial resolution of the DENV case data for each country, we calculated and validated the relevant summary statistics of Index P at the municipality level for Brazil, the province level for Thailand and the country level for Mexico and Puerto Rico rather than at the spatial pixel level. Here, we included several analyses that support the application of Index P to explain past spatiotemporal dynamics of DENV, reaffirming the results of previous studies 5,24,26-32 . Local transmission intensity. To quantify how well the values of Index P characterize spatial variation in local DENV incidence, we compared the mean Index P during a typical year and the mean yearly incidence in each municipality across Brazil and each province across Thailand. A typical yearly time series is estimated by calculating the mean value of a chosen variable (e.g. Index P, log(incidence + 1), incidence, etc.) independently for each of the twelve months of the year across all available years. The map of transmission potential was consistent with the spatial distribution of incidence for both Brazil and Thailand (Fig. 3a,b). In Brazil, the hyperendemic regions in the Northeast and the Midwest had higher estimated transmission potential, while the low incidence areas in the South stretching from Rio Grande do Sul to Minas Gerais had lower estimated transmission potential. Across Thailand the transmission potential was estimated to be uniformly high in line with the relatively high incidence and endemicity of DENV in the country.
Estimated transmission potential and incidence exhibited a nonlinear relationship. Incidence increased with Index P until reaching a plateau in which increases in transmission potential did not yield higher incidence (Fig. 3c). This nonlinear relationship is consistent with the interpretation of Index P as the absolute transmission potential of an adult mosquito, where increases in Index P are tied to increases in R 0 and the probability and scale of epidemic growth 51 . The observed plateau at high Index P values can be explained by the dynamics of viral transmission, where the size of epidemics is constrained by factors other than transmission potential, www.nature.com/scientificdata www.nature.com/scientificdata/ particularly the development of herd immunity over time. As such, areas estimated to have very high transmission potential presented similar levels of high incidence.
The reported incidence was lower than expected from the estimated transmission potential in several municipalities in the Northwest states of Amazonas and Pará in Brazil (Fig. 3a). Given that these areas form part of the Amazon rainforest, these inconsistencies might be due to a combination of poor surveillance, lack of connectivity to allow frequent viral introduction and/or low population density. Indeed, the municipalities within Brazil's Legal Amazon that had lower than expected incidence had some of the largest areas of all municipalities, and relatively low population densities (Fig. 3a,c). Although reported incidence was relatively high across provinces in Thailand, there was some variation with higher than expected incidence in the northwestern provinces (Fig. 3b). This high incidence is reported across recent studies 52-55 and although not entirely understood it is potentially linked to a combination of several factors including variation in reporting rates 55 , influence of long distance 56 and Myanmar cross-border human mobility 52,53 , socio-demographic factors 57,58 and misdiagnosis of other endemic diseases with similar symptomatic profiles 59 .
typical year seasonality. We compared the Index P dynamics to the observed DENV seasonality by calculating Spearman's correlation coefficient between the time series of Index P and log-transformed incidence during a typical year in each municipality in Brazil and each province in Thailand (Fig. 4a,b). Municipalities in Brazil with an average of fewer than 12 cases per year (2725/5570) were excluded since their incidence dynamics did not have sufficient information for a measurable seasonal signal. To quantify the time delay between the Index P and incidence dynamics, we calculated the correlation coefficient for all possible monthly lags (i.e. −5 to 6 months inclusive) and selected the lag with the highest coefficient. The lag with the highest correlation coefficient and the coefficient itself are referred to as the lag and lag-adjusted correlation of a municipality respectively.
The distribution of lag-unadjusted correlation coefficients was left skewed with 42% of municipalities in Brazil and 91% of provinces in Thailand showing relatively high coefficients greater than 0.6 ( Fig. 4a,b). Given the development time of Aedes mosquitoes (~2 weeks), the generation time of DENV in human populations (~2-3 weeks) 60,61 and the time for mosquitoes to respond to climatic signals, it can take time for increases in transmission potential to translate into increases in incidence 5,62,63 . Indeed, the Index P dynamics in 96% of municipalities in Brazil preceded the incidence dynamics by 0-3 months, with the proportion of municipalities with time delays of 0, 1, 2 and 3 months estimated at 8%, 25%, 47% and 15%, respectively (Fig. 4a). Similarly, the Index P dynamics in 95% of provinces in Thailand preceded the incidence dynamics by 0-3 months with the proportion of provinces with time delays of 0, 1, 2 and 3 months estimated at 52%, 32%, 9% and 1%, respectively (Fig. 4b). The mechanisms behind the different distributions of time delays in Brazilian municipalities and Thai Fig. 3 Comparison of reported DENV incidence and mean annual Index P for a typical year. (a,b) Mean yearly incidence in each municipality (left panel) and mean annual Index P for a typical year (right panel) across municipalities in Brazil and provinces in Thailand, respectively. Mean yearly incidence is normalized by the maximum value. (c) The relationship between mean annual Index P and mean yearly incidence (cases per 100,000 population). The area below the lower fence (Q1-1.5 × IQR for (0, 0.5) and (0.5, ∞)) which encompasses outliers is shaded grey with municipalities that form part of Brazil's Legal Amazon highlighted black (these municipalities are also highlighted in panel a).
www.nature.com/scientificdata www.nature.com/scientificdata/ provinces could not be identified, but such differences possibly reflect a combination of multiple factors such as local variation in mosquito populations or accumulated herd immunity. Of particular note is that Thailand has historically been endemic to dengue, while Brazil has only recently become endemic to all four serotypes. When we corrected for this time delay, the dynamics for Index P and incidence were highly correlated with 82% of the municipalities in Brazil and 88% of the provinces in Thailand with an average of 12 or more cases per year showing correlation coefficients greater than 0.90 (Fig. 4a,b).
Some municipalities in Brazil (4%) had extreme lags (i.e. negative or greater than 3 months) inconsistent with the biological interpretation of Index P (in the North, spanning the states of Roraima, Pará, Amapá and Maranhão, Fig. 4a). Given that these regions have monsoon-driven climates characterized by intense timely annual rainfall, a factor not considered in Index P but known to influence mosquito abundance, precipitation may play an outside role in determining seasonality in these regions. www.nature.com/scientificdata www.nature.com/scientificdata/ Long term seasonality. We also compared the time series of Index P to the reported DENV cases over many years at the country level for Brazil, Mexico, Puerto Rico and Thailand (Fig. 4c). The lag-unadjusted correlation coefficients were 0.32, 0.37, 0.34 and 0.58 for Brazil, Mexico, Puerto Rico and Thailand, respectively. Consistent with the analysis of typical year seasonality, the Index P dynamics preceded the observed DENV seasonality by 1 or 2 months. Once we adjusted for the predicted time delay, the correlation coefficients increased to 0.72, 0.62, 0.47 and 0.61, respectively.

Usage Notes
Researchers can use Index P maps to explore monthly transmission potential of DENV spread by Ae. aegypti mosquitoes at the spatial pixel level (~28 km 2 ) from 1981 to 2019. There are, however, several caveats to the application of these maps. Firstly, from our technical validation, it can be concluded that areas with low estimated transmission potential can be expected to have lower incidence. However, areas with high estimated potential are susceptible to more uncertainty as non-climatic factors, particularly accumulated herd immunity, human/vector densities and human behavior, can influence whether and how favorable climatic conditions translate into higher incidence. High Index P thus should not be interpreted as a guarantee of transmission but rather as one indicator of favorable climatic conditions for DENV circulation. Users are encouraged to interpret the magnitude of Index P in a region of interest within the context of the social, demographic and ecological conditions that may influence the interaction between climate and local dengue transmission. Secondly, while Index P dynamics reliably capture the seasonality of dengue transmission, there is likely to be a time delay between the proposed climatic effects and the observed changes in incidence, the length of which can vary between locations. Our analyses, and those of others 62,63 , suggest that this time delay is between about 0 to 3 months, but a comparison of past incidence and Index P dynamics at a local level is necessary for more precise, locally-specific estimates. Finally, Index P is unlikely to mimic incidence patterns at very high spatiotemporal resolutions due to obfuscation of the seasonal signals of incidence by factors such as imperfect surveillance and local stochasticity. Thus Index P is most informative at aggregated spatiotemporal resolutions where incidence is expected to demonstrate clear seasonal signals.
To encourage reuse and deepen understanding of the Index P methodology, we developed a new version of the Mosquito-borne Viral Suitability Estimator (MVSE) software package for the R programming environment 24 . This R package provides a set of related functions that can be used to estimate Index P time series given temperature and humidity time series and user-defined probability distributions of the biological parameters for the selected host/vector/virus system. A short tutorial on the features of the MVSE package is also made available. Furthermore, we have provided an R Markdown document that highlights useful R software packages and functions for the visualization and analysis of the dataset of Index P maps. Though all code used to generate this dataset was written in R and C++, the maps themselves are provided as TIFF files to facilitate cross-compatibility.

code availability
The MVSE R package can be installed from GitHub (https://github.com/TaishiNakase/MVSE). The MVSE package tutorial and the R Markdown document that provides example code for the visualization of the Index P maps can also be found on GitHub (https://github.com/TaishiNakase/Index-P-estimation-and-applications).